function [ volume, surfacearea, areaAbove, areaBelow ] = calcVolumeSurfacearea( contour )
%% barycenter
[ geom, ~, ~ ] = polygeom( contour(:,1), contour(:,2) );
yb = geom(3);

% figure(1);
% plot(contour(:,1),contour(:,2),'.-', [min(contour(:,1)) max(contour(:,1))],[yb yb],'-');

%% distinguish points above and below symmetry line
maskAbove = contour(:,2) >= yb;
idxCircShift = find(~[maskAbove(end); maskAbove(1:end-1)] & maskAbove, 1);
maskAbove = circshift(maskAbove, [-(idxCircShift-1) 0]);
contour = circshift(contour, [-(idxCircShift-1) 0]);

contourAbove = contour(maskAbove,:);
contourBelow = contour(~maskAbove,:);

% figure(1);
% plot(contourAbove(:,1),contourAbove(:,2),'.-', contourBelow(:,1),contourBelow(:,2),'.-', [min(contour(:,1)) max(contour(:,1))],[yb yb],'-');

%% insert intersection points with symmetry axis
pt1 = contourAbove(1,:);
pt2 = contourBelow(end,:);
xIntersect1 = (yb-pt1(2))*(pt2(1)-pt1(1))/(pt2(2)-pt1(2))+pt1(1);
contourAbove = [xIntersect1 yb; contourAbove];
contourBelow = [contourBelow; xIntersect1 yb];

pt1 = contourAbove(end,:);
pt2 = contourBelow(1,:);
xIntersect2 = (yb-pt1(2))*(pt2(1)-pt1(1))/(pt2(2)-pt1(2))+pt1(1);
contourAbove = [contourAbove; xIntersect2 yb];
contourBelow = [xIntersect2 yb; contourBelow];

% figure(1);
% plot(contourAbove(:,1),contourAbove(:,2),'.-', contourBelow(:,1),contourBelow(:,2),'.-', [min(contour(:,1)) max(contour(:,1))],[yb yb],'-');

%% calculate volume and surface area
[ geom, ~, ~ ] = polygeom( contourAbove(:,1), contourAbove(:,2) );
areaAbove = geom(1);
radiusAbove = abs(geom(3) - yb);
VolumeAbove = areaAbove * 2*pi*radiusAbove;

[lengthAbove, ~, posAbove] = calcOpenPolygon( contourAbove(:,1), contourAbove(:,2) );
radiusAbove = abs(posAbove - yb);
SurfaceareaAbove = lengthAbove * 2*pi*radiusAbove;

[ geom, ~, ~ ] = polygeom( contourBelow(:,1), contourBelow(:,2) );
areaBelow = geom(1);
radiusBelow = abs(geom(3) - yb);
VolumeBelow = areaBelow * 2*pi*radiusBelow;

[lengthBelow, ~, posBelow] = calcOpenPolygon( contourBelow(:,1), contourBelow(:,2) );
radiusBelow = abs(posBelow - yb);
SurfaceareaBelow = lengthBelow * 2*pi*radiusBelow;

volume = 0.5*(abs(VolumeAbove) + abs(VolumeBelow));
surfacearea = 0.5*(SurfaceareaAbove + SurfaceareaBelow);
end